subroutine solicitar(it,tolerancia,a,b)
	implicit none
	integer(kind=4)::it
	real(kind=8)::tolerancia,a,b
	print*,'Inserte el numero de iteraciones.'
	read*,it
	print*,'Inserte la tolerancia de error.'
	read'(e7.5)',tolerancia
	print*,'Inserte los dos valores del intervalo.'	
	read*,a,b
end subroutine solicitar

function funcion(x,A1,A2) result(x2)
	implicit none
	real(kind=8),intent(in)::x,A1,A2
   real(kind=8)::x2
	x2=A1-x*A2
end function funcion

function leer(vectx,vectw) result(n)
	implicit none
	real(kind=8),dimension(:),allocatable::vectx,vectw
	integer(kind=4)::n
	open(unit=1,file='fichero.dat',status='old')
	rewind(1)  		
	read(1,*)n
	!print*,n
	allocate(vectx(0:n))
	allocate(vectw(0:n))
	read(1,*) vectx
	read(1,*) vectw
   !print*,vectx
   !print*,vectw
	close(1)
end function leer

function leerAplicacion(vectx,vectw,ea) result(n)
	implicit none
	real(kind=8),dimension(:),allocatable::vectx,vectw
	integer(kind=4)::n
	real(kind=8)::ea
	open(unit=1,file='fichero.dat',status='old')
	rewind(1)  		
	read(1,*)ea
   read(1,*)n 	
	allocate(vectx(n))
	allocate(vectw(n))
	read(1,*) vectx
	read(1,*) vectw
	close(1)
end function leerAplicacion



function solicitarPtos(ptos) result(p)
	implicit none
	integer(kind=4)::p,i
	real(kind=8),dimension(:),allocatable::ptos
	print*,'¿Cuántos puntos desea evaluar?'
	read*,p
	allocate(ptos(p))
	print*,'Inserte los puntos de evaluacion: '	
	do i=1,p
		read*,ptos(i)
	end do
end function solicitarPtos

function difDiv(vectx,vectw,n) result(coef)
	implicit none
	integer(kind=4)::n,i,j,k
	real(kind=8),dimension(:),allocatable::vectx,vectw,coef
	real(kind=8),dimension(:,:),allocatable::tabla
	allocate(tabla(0:n,0:n))
	!do i=1,n
  	!	do j=n,i,-1
    !		vectw(j)=(vectw(j)-vectw(j-1))/(vectx(j)-vectx(j-i))
  		!end do
!	end do
!	print*,vectw	
	allocate(coef(0:n))
	open(unit=2,file='DD.dat',status='unknown')
	tabla=0
	do i=0,n
		tabla(i,0)=vectw(i)
	end do
	do j=1,n
		do i=0,n-j
			tabla(i,j)=(tabla(i,j-1)-tabla(i+1,j-1))/(vectx(i)-vectx(i+j))	
		end do
	end do
	do i=0,n
		write(2,*) tabla(i,:)	
	end do
	coef=tabla(0,0:n)
	close(2)
end function difDiv

function lagrange(n,vectx,coef,x) result(x2)
	integer(kind=4)::n,i,j
	real(kind=8),dimension(:),allocatable::vectx,coef
	real(kind=8)::x,x2
	x2=coef(n)
	do i=n-1,0,-1
		x2=coef(i)+(x-vectx(i))*x2
	end do 
end function lagrange

function eleccion() result(caso)
	integer(kind=4)::caso,cuad,func
	print*,'Desea usar Trapecio o Simpson (1/2)'
	read*,cuad
	print*,'Desde fichero o funcion conocida(1/2)'
	read*,func
	if (cuad==1) then
		if (func==1) then
			caso=1
		elseif (func==2) then 
			caso=2
		end if
	elseif (cuad==2) then
		if (func==1) then
			caso=3
		elseif (func==2) then
			caso=4
		end if
	end if
end function eleccion


subroutine dicotomia(it,tolerancia,a,b,A1,A2)
   implicit none   
   integer(kind=4)::it,i
   real(kind=8)::tolerancia,a,b,x,x0,fx,fa,A1,A2
   interface   	
	   function funcion(x,A1,A2) result(x2)
         real(kind=8),intent(in)::x,A1,A2
         real(kind=8)::x2
      end function funcion
   end interface
	x=0._8
	if (funcion(a,A1,A2)==0) then
		x=a
	else if (funcion(b,A1,A2)==0) then
		x=b
	end if	
	if (funcion(a,A1,A2)*funcion(b,A1,A2)<0) then
		print*,'entro3'
   	open(unit=1,file='dicotomia.dat',status='unknown')
   	write(1,'(a14,i3,a14,e12.6)') "iteraciones:  ",it,"  tolerancia:  ",tolerancia
		write(1,'(a4,f19.14,a4,f19.14)') "a: ",a," b: ",b
		x=a
 	   do i=1,it	
			x0=x
			x=(a+b)*0.5_8
			fx=funcion(x,A1,A2)
			fa=funcion(a,A1,A2)
			if (fx==0) then
				exit
			end if			
			if((fa*fx)<0) then
				b=x
			else
				a=x
			end if
   	   write(1,'(i3,f19.14,e15.7)') i,x,abs(x-x0)/abs(x)
			if (abs(x-x0)<tolerancia*abs(x)) then
   	      exit
   	   end if  
  	  end do
  	  close(1)
   end if
	print*,x
end subroutine dicotomia

function s_trapecio(a,b,fa,fb) result(t)
	implicit none
	real(kind=8),intent(in)::a,b,fa,fb
   real(kind=8)::t
	t=(b-a)*(fa+fb)*0.5_8
!	print*,t
	
end function s_trapecio

function s_simpson(a,b,fa,fb,fc) result(t)
	implicit none
	real(kind=8),intent(in)::a,b,fa,fb,fc
   real(kind=8)::t
	t=((b-a)/6._8)*(fa+4._8*fc+fb)
!	print*,t
end function s_simpson

function Trapecio(a,b,interv,val,div) result(t)
	integer(kind=4)::interv,i			
	real(kind=8)::a,b,x,div
	real(kind=8),dimension(:),allocatable::val
	real(kind=8)::t
	x=a
	t=0._8
	do i=1,interv
		t=t+s_trapecio(x,x+div,val(i),val(i+1))
		x=x+div
	end do
	print*,t
end function Trapecio

function TrapecioFichero(vectx,vectw,n) result(t)
	integer(kind=4)::n,i			
!	real(kind=8)::a,b,x,div
	real(kind=8),dimension(:),allocatable::vectx,vectw
	real(kind=8)::t
	t=0._8
	do i=0,n-1
		t=t+s_trapecio(vectx(i),vectx(i+1),vectw(i),vectw(i+1))
	end do
	print*,t
end function TrapecioFichero

function Simpson(a,b,interv,val,div) result(t)
	integer(kind=4)::interv			
	real(kind=8)::a,b,x,div
	real(kind=8),dimension(:),allocatable::val
	real(kind=8)::t
	t=0._8
	x=a
	do i=1,2*interv,2
		t=t+s_simpson(x,x+2*div,val(i),val(i+2),val(i+1))
		x=x+2*div	
	end do
	print*,t
end function Simpson
	
function SimpsonFichero(vectx,vectw,vectmedio,n) result(t)
	integer(kind=4)::i,n		
	real(kind=8),dimension(:),allocatable::vectx,vectw,vectmedio
	real(kind=8)::t
	t=0._8
	do i=0,n-1
		t=t+s_simpson(vectx(i),vectx(i+1),vectw(i),vectw(i+1),vectmedio(i))
	end do
	print*,t
end function SimpsonFichero

program practica3
	implicit none	
	integer(kind=4)::it,n,p,i,caso
	real(kind=8)::tolerancia,a,b,c,integral,eval,div,ea,A1,A2
	real(kind=8),dimension(:),allocatable::vectx,vectw,ptos,coef,vectmedio,vectcos,vectfunc
	interface	
	   function funcion(x,A1,A2) result(x2)
         real(kind=8),intent(in)::x,A1,A2
         real(kind=8)::x2
      end function funcion
		subroutine dicotomia(it,tolerancia,a,b,A1,A2)
			integer(kind=4)::it			
			real(kind=8)::tolerancia,a,b,A1,A2
		end subroutine dicotomia
		subroutine solicitar(it,tolerancia,a,b)
			integer(kind=4)::it
			real(kind=8)::tolerancia,a,b
		end subroutine solicitar
		function s_trapecio(a,b,fa,fb) result(t)
			real(kind=8),intent(in)::a,b,fa,fb
   		real(kind=8)::t
		end function s_trapecio
		function s_simpson(a,b,fa,fb,fc) result(t)
			real(kind=8),intent(in)::a,b,fa,fb,fc
   		real(kind=8)::t
		end function s_simpson
		function leer(vectx,vectw) result(n)
			real(kind=8),dimension(:),allocatable::vectx,vectw
			integer(kind=4)::n
		end function leer
		function leerAplicacion(vectx,vectw,ea) result(n)
			implicit none
			real(kind=8),dimension(:),allocatable::vectx,vectw
			integer(kind=4)::n
			real(kind=8)::ea
		end function leerAplicacion
		function solicitarPtos(ptos) result(p)
			integer(kind=4)::p
			real(kind=8),dimension(:),allocatable::ptos
		end function solicitarPtos
		function difDiv(vectx,vectw,n) result(coef) 
			integer(kind=4)::n
			real(kind=8),dimension(:),allocatable::vectx,vectw,coef
		end function difDiv
		function lagrange(n,vectx,coef,x) result(x2)
			integer(kind=4)::n
			real(kind=8),dimension(:),allocatable::vectx,coef
			real(kind=8)::x,x2				
		end function lagrange
		function eleccion() result(caso)
			integer(kind=4)::caso
		end function eleccion
		function Trapecio(a,b,interv,val,div) result(t)
			integer(kind=4)::interv			
			real(kind=8)::a,b,div
			real(kind=8),dimension(:),allocatable::val
   		real(kind=8)::t
		end function Trapecio
		function Simpson(a,b,interv,val,div) result(t)
			integer(kind=4)::interv			
			real(kind=8)::a,b,div
			real(kind=8),dimension(:),allocatable::val
   		real(kind=8)::t
		end function Simpson		
		function TrapecioFichero(vectx,vectw,n) result(t)
			integer(kind=4)::n			
			real(kind=8),dimension(:),allocatable::vectx,vectw
			real(kind=8)::t
		end function TrapecioFichero
		function SimpsonFichero(vectx,vectw,vectmedio,n) result(t)
			integer(kind=4)::n		
			real(kind=8),dimension(:),allocatable::vectx,vectw,vectmedio
			real(kind=8)::t
		end function SimpsonFichero
	end interface
	!call solicitar(it,tolerancia,a,b)
	!call dicotomia(it,tolerancia,a,b)
	!integral=s_trapecio(0._8,2._8,funcion(0._8),funcion(2._8))
	!integral=s_simpson(-14._8,16._8,funcion(-14._8),funcion(16._8),funcion((-14._8+16._8)*0.5_8))
	!print*,integral
	!n=leer(vectx,vectw)
	!allocate(coef(0:n))
	!coef=difDiv(vectx,vectw,n)
	!p=solicitarPtos(ptos)
	!do i=1,p
	!	eval=lagrange(n,vectx,coef,ptos(i))
	!	print*,eval
	!end do
	
	
	n=leerAplicacion(vectx,vectw,ea)
	allocate(vectfunc(n))
	allocate(vectcos(n))
	allocate(coef(0:n))
	coef=difDiv(vectx,vectw,n)	
	do i=1,n
		vectfunc(i)=500000*ea*vectw(i)*3.141592/180._8
		vectcos(i)=30*cos(vectx(i)*(3.141592/180._8))		
	end do
	!print*,vectfunc
	!print*,vectcos
	A1=TrapecioFichero(vectx,vectfunc,n)
	A2=TrapecioFichero(vectx,vectcos,n)
	call solicitar(it,tolerancia,a,b)
	call dicotomia(it,tolerancia,vectx(1),vectx(n),A1,A2)	
	print*,funcion(5.24469506149841_8,A1,A2)
!	caso=eleccion()
	select case(caso)
		case(1) !leer fichero y usar Trapecio
			n=leer(vectx,vectw)
			!eval=Trapecio(vectx(0),vectx(n-1),n,vectw,div)
			eval=TrapecioFichero(vectx,vectw,n)
		case(2) !usar funcion y Trapecio
			print*,'Introduzca los extremos del intervalo'
			read*,a,b
			print*,'Introduzca el numero de divisiones'
			read*,n
			div=(b-a)/n!OJO con esta n y con vectw
			allocate(vectw(n+1))
			c=a
			do i=1,n+1
				vectw(i)=funcion(c,1._8,2._8)
				c=c+div
			end do
			eval=Trapecio(a,b,n,vectw,div)
		case(3) !leer fichero y usar Simpson
			n=leer(vectx,vectw)
			allocate(coef(0:n))
			coef=difDiv(vectx,vectw,n)
			allocate(vectmedio(0:n-1))			
			do i=0,n-1
				vectmedio(i)=lagrange(n,vectx,coef,(vectx(i+1)+vectx(i))/2)
			end do
			eval=SimpsonFichero(vectx,vectw,vectmedio,n)

		case(4) !usar funcion y Simpson
			print*,'Introduzca los extremos del intervalo'
			read*,a,b
			print*,'Introduzca el numero de divisiones'
			read*,n
			div=(b-a)/n
			div=div/2
			allocate(vectw(2*n+1))
			c=a
			do i=1,2*n+1
				vectw(i)=funcion(c,1._8,2._8)
				c=c+div
			end do
			eval=Simpson(a,b,n,vectw,div)
	end select
end program practica3
